#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")
In this step, we are loading the libraries which are required for analysis and plotting of the data.
library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(limma)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)
my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")
# Mouse 1
saline_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/Saline_sample_round_one/outs/"
list.files(saline_mouse_1)
## [1] "analysis" "cloupe.cloupe"
## [3] "filtered_feature_bc_matrix" "filtered_feature_bc_matrix.h5"
## [5] "metrics_summary.csv" "molecule_info.h5"
## [7] "possorted_genome_bam.bam" "possorted_genome_bam.bam.bai"
## [9] "probe_set.csv" "raw_feature_bc_matrix"
## [11] "raw_feature_bc_matrix.h5" "spatial"
## [13] "spatial_enrichment.csv" "web_summary.html"
saline_mouse_1_data <- Load10X_Spatial(
saline_mouse_1,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "saline",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 2
saline_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/Saline_sample_round_two/outs/"
list.files(saline_mouse_2)
## [1] "analysis" "cloupe.cloupe"
## [3] "filtered_feature_bc_matrix" "filtered_feature_bc_matrix.h5"
## [5] "metrics_summary.csv" "molecule_info.h5"
## [7] "possorted_genome_bam.bam" "possorted_genome_bam.bam.bai"
## [9] "probe_set.csv" "raw_feature_bc_matrix"
## [11] "raw_feature_bc_matrix.h5" "saline_web_summary.html"
## [13] "spatial" "spatial_enrichment.csv"
saline_mouse_2_data <- Load10X_Spatial(
saline_mouse_2,
filename = "filtered_feature_bc_matrix.h5",
assay = "spatial",
slice = "saline",
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL)
# Mouse 1
saline_mouse_1_data <- SCTransform(saline_mouse_1_data, assay = "spatial", verbose = FALSE)
# Mouse 2
saline_mouse_2_data <- SCTransform(saline_mouse_2_data, assay = "spatial", verbose = FALSE)
p1 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
p1 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
# Mouse 1
saline_mouse_1_analyzed <- RunPCA(saline_mouse_1_data, assay = "SCT", verbose = FALSE)
saline_mouse_1_analyzed <- FindNeighbors(saline_mouse_1_analyzed, reduction = "pca", dims = 1:5)
saline_mouse_1_analyzed <- FindClusters(saline_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
saline_mouse_1_analyzed <- RunUMAP(saline_mouse_1_analyzed, reduction = "pca", dims = 1:5)
# Mouse 2
saline_mouse_2_analyzed <- RunPCA(saline_mouse_2_data, assay = "SCT", verbose = FALSE)
saline_mouse_2_analyzed <- FindNeighbors(saline_mouse_2_analyzed, reduction = "pca", dims = 1:5)
saline_mouse_2_analyzed <- FindClusters(saline_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
saline_mouse_2_analyzed <- RunUMAP(saline_mouse_2_analyzed, reduction = "pca", dims = 1:5)
# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)
# Mouse 1
saline_mouse_1_analyzed.sce <- as.SingleCellExperiment(saline_mouse_1_analyzed)
# Mouse 2
saline_mouse_2_analyzed.sce <- as.SingleCellExperiment(saline_mouse_2_analyzed)
#create reference data
ref.data <- celldex::ImmGenData()
# run singleR pipeline to find the cell types
## Mouse 1
cell_types_mouse_1 <- SingleR(test=saline_mouse_1_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
## Mouse 2
cell_types_mouse_2 <- SingleR(test=saline_mouse_2_analyzed.sce, assay.type.test=1,
ref=ref.data, labels=ref.data$label.main)
# view cell types
# Mouse 1
table(cell_types_mouse_1$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 58 41 1095 262
## Fibroblasts ILC Macrophages Monocytes
## 384 1 1007 19
## Neutrophils Stromal cells
## 2 180
#checking cell types Mouse 1
plotScoreHeatmap(cell_types_mouse_1)
plotDeltaDistribution(cell_types_mouse_1)
summary(is.na(cell_types_mouse_1$pruned.labels))
## Mode FALSE
## logical 3049
# Mouse 2
table(cell_types_mouse_2$labels)
##
## B cells DC Endothelial cells Epithelial cells
## 60 11 706 62
## Fibroblasts Macrophages Monocytes Stromal cells
## 162 1534 9 35
#checking cell types Mouse 2
plotScoreHeatmap(cell_types_mouse_2)
plotDeltaDistribution(cell_types_mouse_2)
summary(is.na(cell_types_mouse_2$pruned.labels))
## Mode FALSE TRUE
## logical 2570 9
# Mouse 1
saline_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels
# Mouse 2
saline_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels
p1 <- DimPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
p1 <- DimPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)
p2 <- SpatialDimPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)
p1 + p2
VlnPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
VlnPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols = my_cols)
B_cell_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Ighd 8.255271e-64 1.3507961 0.828 0.118 1.269991e-59
## Cd19 8.357792e-61 1.4494167 0.897 0.155 1.285763e-56
## Il9r 8.629871e-57 0.4171455 0.345 0.018 1.327619e-52
## Ms4a1 2.148405e-52 0.7874505 0.534 0.053 3.305106e-48
## Bank1 3.460125e-52 0.7886827 0.603 0.069 5.323056e-48
## Cr2 3.013761e-50 0.8212669 0.500 0.048 4.636370e-46
## Tnfrsf13c 6.413623e-45 1.0369011 0.707 0.117 9.866717e-41
## Blnk 7.485194e-45 0.9015819 0.741 0.126 1.151522e-40
## Spib 1.319907e-44 1.1623285 0.707 0.118 2.030545e-40
## Fcrla 1.453135e-44 0.6200531 0.483 0.050 2.235503e-40
write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/B_cell_mouse_1.csv")
DC_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Zbtb46 2.592652e-08 0.3017862 0.366 0.098 0.0003988535
## Gbp4 3.439860e-08 0.5729455 0.927 0.584 0.0005291880
## Epas1 6.949482e-08 -0.8047154 0.512 0.870 0.0010691084
## Fmo2 1.405798e-07 -0.8504219 0.341 0.683 0.0021626804
## Dffa 1.742615e-07 0.3038923 0.415 0.129 0.0026808388
## Inmt 5.366006e-07 -0.8811450 0.488 0.840 0.0082550640
## Nos2 5.879765e-07 0.5272287 0.902 0.527 0.0090454306
## Rogdi 3.931562e-06 0.3697050 0.634 0.295 0.0604831477
## Lyrm2 4.272549e-06 0.2744930 0.317 0.098 0.0657289010
## Man1c1 4.734111e-06 0.2995503 0.415 0.151 0.0728295623
write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/DC_mouse_1.csv")
Stromal_cells_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igkc 2.148446e-18 -0.5779135 1.000 1.000 3.305170e-14
## Igha 2.454418e-18 -0.5396479 1.000 0.998 3.775877e-14
## Gsn 2.048482e-15 0.5022576 0.911 0.738 3.151385e-11
## Inmt 8.253679e-14 0.5389418 0.961 0.827 1.269746e-09
## Ly6e 1.938211e-13 -0.3734319 0.956 0.985 2.981744e-09
## Tnfaip2 5.132503e-13 -0.4556705 0.917 0.966 7.895842e-09
## Jchain 1.233062e-12 -0.5806728 0.650 0.864 1.896942e-08
## Ankrd1 3.898149e-10 0.5452627 0.383 0.199 5.996912e-06
## H2-DMb2 1.064259e-09 -0.3793610 0.883 0.931 1.637256e-05
## Ighj1 1.381416e-09 -0.3043841 0.656 0.845 2.125171e-05
write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Stromal_cells_mouse_1.csv")
Endothelial_cells_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cldn5 1.027986e-53 0.5591739 0.916 0.772 1.581453e-49
## Igfbp2 3.888415e-48 0.6519764 0.915 0.801 5.981938e-44
## Epas1 8.744963e-45 0.4787279 0.934 0.827 1.345325e-40
## Rasip1 2.717011e-40 0.3987670 0.694 0.464 4.179849e-36
## Tns1 1.906312e-39 0.4303160 0.953 0.872 2.932671e-35
## Timp3 4.363481e-37 0.3786639 0.988 0.942 6.712779e-33
## Clec14a 4.322131e-36 0.3560302 0.554 0.328 6.649166e-32
## Kdr 8.719974e-34 0.3041826 0.526 0.305 1.341481e-29
## Tmcc2 1.355709e-33 0.3613636 0.595 0.375 2.085622e-29
## Igkc 6.254199e-33 -0.3852192 1.000 0.999 9.621460e-29
write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Endothelial_cells_mouse_1.csv")
Epithelial_cell_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Inmt 1.143840e-19 -0.6655945 0.706 0.847 1.759683e-15
## Cldn5 5.345920e-19 -0.5996306 0.691 0.836 8.224164e-15
## Tns1 4.602813e-16 -0.4876424 0.821 0.908 7.080967e-12
## Ltbp4 5.284907e-16 -0.4639505 0.206 0.480 8.130301e-12
## Timp3 5.455946e-16 -0.4549849 0.905 0.964 8.393427e-12
## Igha 1.112342e-15 0.3166163 1.000 0.998 1.711226e-11
## Igkc 3.078809e-15 0.3298940 1.000 1.000 4.736440e-11
## Pecam1 6.891285e-14 -0.3491782 0.179 0.435 1.060155e-09
## Ramp2 1.604636e-13 -0.3287006 0.092 0.312 2.468573e-09
## Ager 2.040610e-13 -0.3179175 0.973 0.991 3.139274e-09
write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm+bcg/Epithelial_cell_mouse_1.csv")
Fibroblasts_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Inmt 5.266062e-22 0.4666610 0.966 0.816 8.101309e-18
## Gsn 1.647372e-17 0.4518127 0.878 0.730 2.534317e-13
## Nbl1 1.612737e-14 0.3469147 0.859 0.757 2.481034e-10
## Clec3b 9.523951e-14 0.4607619 0.576 0.420 1.465165e-09
## Fmo2 9.902154e-11 0.2931271 0.794 0.661 1.523347e-06
## Ctss 1.366486e-10 -0.2699602 0.984 0.991 2.102203e-06
## Ltbp4 2.049154e-10 0.3284453 0.581 0.438 3.152419e-06
## Ptprc 1.829641e-09 -0.2610686 0.646 0.768 2.814719e-05
## Rarres2 1.653134e-08 0.3535110 0.729 0.629 2.543181e-04
## Col1a2 1.811325e-08 0.2753930 0.831 0.737 2.786543e-04
write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Fibroblasts_mouse_1.csv")
Macrophages_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_1, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Igha 1.802204e-45 0.3570873 1.000 0.997 2.772510e-41
## Inmt 2.349932e-42 -0.5859870 0.761 0.872 3.615135e-38
## Igfbp2 2.445557e-42 -0.6746716 0.776 0.875 3.762244e-38
## Igkc 3.043211e-37 0.2814937 1.000 1.000 4.681676e-33
## Epas1 4.994383e-35 -0.4661141 0.805 0.895 7.683358e-31
## Cldn5 6.077221e-35 -0.5096620 0.761 0.855 9.349197e-31
## Timp3 6.586934e-33 -0.4033961 0.944 0.966 1.013334e-28
## Tns1 4.589537e-32 -0.4275692 0.865 0.918 7.060544e-28
## Ctss 1.098293e-31 0.3599147 0.998 0.987 1.689614e-27
## Tnfaip2 1.227430e-31 0.3164347 0.995 0.947 1.888278e-27
write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophages_mouse_1.csv")
B_cell_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cxcr5 1.397901e-79 1.2684368 0.817 0.089 1.908834e-75
## Bank1 2.629509e-75 1.6224021 0.950 0.146 3.590595e-71
## Cr2 4.492997e-68 1.4266862 0.700 0.077 6.135188e-64
## Il9r 8.464439e-64 0.6025927 0.467 0.031 1.155819e-59
## Pax5 3.205187e-63 1.1707456 0.733 0.090 4.376683e-59
## Fcrl5 5.607384e-63 0.9722059 0.600 0.057 7.656883e-59
## Dtx1 3.669356e-61 0.6956637 0.550 0.047 5.010505e-57
## Ccr6 8.135732e-61 1.9387029 0.950 0.200 1.110934e-56
## Tnfrsf13c 1.292829e-58 1.8747025 0.917 0.178 1.765358e-54
## Ighd 6.070972e-58 1.8733801 0.883 0.167 8.289912e-54
write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/B_cell_mouse_2.csv")
DC_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Klhdc9 1.052979e-23 0.3384042 0.273 0.007 1.437843e-19
## Clec9a 5.968194e-19 0.4251603 0.364 0.016 8.149569e-15
## Cfap74 8.906069e-19 0.5914211 0.455 0.025 1.216124e-14
## Mst1 1.538070e-14 0.3311671 0.273 0.012 2.100234e-10
## Zfp160 1.356914e-13 0.4892451 0.455 0.035 1.852866e-09
## Erich3 1.016938e-11 1.0236002 0.545 0.059 1.388628e-07
## Lrrc43 4.023684e-11 0.9656281 0.636 0.086 5.494340e-07
## March4 5.542413e-11 0.4979459 0.364 0.028 7.568166e-07
## Fzd8 1.651961e-10 0.3223094 0.273 0.017 2.255753e-06
## Gp2 6.170181e-10 0.4113967 0.273 0.018 8.425383e-06
write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/DC_mouse_2.csv")
Stromal_cells_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Myh11 2.799853e-12 0.9352038 0.571 0.163 3.823199e-08
## Cst3 4.295150e-12 0.6681956 1.000 1.000 5.865027e-08
## Des 7.746254e-12 1.0718122 0.857 0.488 1.057751e-07
## Tagln 7.373254e-11 0.9235898 0.771 0.317 1.006818e-06
## Nbl1 4.602124e-10 0.7027256 1.000 0.985 6.284200e-06
## Nphp4 5.992370e-10 0.3655504 0.314 0.060 8.182581e-06
## Gsn 8.847970e-10 0.8142324 1.000 0.927 1.208190e-05
## Fth1 8.466740e-09 -0.4846324 1.000 1.000 1.156133e-04
## Pcp4l1 9.268633e-09 0.4812146 0.457 0.126 1.265632e-04
## Col14a1 1.227715e-08 0.3198091 0.314 0.066 1.676445e-04
write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Stromal_cells_mouse_2.csv")
Endothelial_cells_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Cldn5 1.411386e-98 0.8332450 0.986 0.896 1.927248e-94
## Timp3 2.632595e-76 0.5422023 1.000 0.996 3.594809e-72
## Epas1 1.200408e-75 0.6283701 0.992 0.953 1.639157e-71
## Acer2 2.821727e-67 0.5590838 0.993 0.952 3.853068e-63
## Igfbp2 4.746039e-67 1.0716958 0.877 0.690 6.480717e-63
## Egfl7 1.783598e-66 0.5758111 0.992 0.915 2.435503e-62
## Ace 4.473554e-62 0.6363960 0.967 0.874 6.108638e-58
## Tspan7 6.054914e-62 0.6679605 0.856 0.612 8.267984e-58
## Plvap 1.008573e-58 0.4802857 0.997 0.981 1.377207e-54
## Pecam1 3.042090e-56 0.5633490 0.939 0.798 4.153974e-52
write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Endothelial_cells_mouse_2.csv")
Epithelial_cell_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1700016K19Rik 4.489550e-59 0.9366679 0.629 0.068 6.130480e-55
## Dynlrb2 5.376514e-48 0.7171786 0.484 0.048 7.341630e-44
## 1700013F07Rik 1.266140e-47 0.5067919 0.387 0.030 1.728915e-43
## Wdr63 1.681633e-45 0.6232098 0.484 0.050 2.296271e-41
## Kndc1 4.331139e-44 0.6357987 0.468 0.049 5.914171e-40
## March4 7.567509e-43 0.4041338 0.323 0.023 1.033343e-38
## Tekt1 2.856782e-41 1.0867857 0.629 0.101 3.900935e-37
## Rp1 3.613835e-41 1.0314188 0.629 0.100 4.934692e-37
## Alox12e 5.386355e-41 0.3617788 0.274 0.017 7.355067e-37
## Dnah5 1.044865e-39 0.7036400 0.532 0.072 1.426763e-35
write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Epithelial_cell_mouse_2.csv")
Fibroblasts_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Col1a2 3.272175e-17 0.4468853 1.000 0.998 4.468154e-13
## Igha 5.584053e-15 1.5201650 0.994 1.000 7.625025e-11
## C3 1.046965e-14 0.2901238 1.000 1.000 1.429631e-10
## Apoe 6.533885e-14 0.4792034 0.957 0.892 8.922020e-10
## Gbp2b 2.702983e-13 -0.3358854 1.000 1.000 3.690923e-09
## Igkc 6.114446e-13 1.3331430 0.994 1.000 8.349276e-09
## Col1a1 6.408404e-13 0.4230387 0.988 0.969 8.750676e-09
## Ctss 9.100994e-13 -0.3472002 1.000 1.000 1.242741e-08
## Mzb1 2.642364e-12 0.7704413 0.691 0.467 3.608148e-08
## Col3a1 5.232022e-12 0.4323003 1.000 0.994 7.144325e-08
write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Fibroblasts_mouse_2.csv")
Macrophages_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Macrophages", min.pct = 0.25)
head(Macrophages_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Ctss 9.097211e-77 0.4498103 1.000 1.000 1.242224e-72
## Fth1 7.823133e-71 0.3583520 1.000 1.000 1.068249e-66
## Psap 5.709096e-56 0.3176438 1.000 1.000 7.795771e-52
## Cfb 1.370117e-55 0.4771078 0.994 0.988 1.870895e-51
## Itgb2 4.442422e-51 0.4104689 0.998 0.984 6.066127e-47
## Il1rn 1.955168e-49 0.5930291 0.984 0.958 2.669781e-45
## Lgals3 3.729628e-48 0.4458466 0.999 0.995 5.092807e-44
## Nos2 1.335862e-47 0.5751555 0.977 0.939 1.824119e-43
## Cyba 2.697403e-46 0.2925800 1.000 1.000 3.683303e-42
## Ctsb 7.022225e-42 0.3247272 1.000 1.000 9.588848e-38
write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophages_mouse_2.csv")
Monocytes_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "Monocytes", min.pct = 0.25)
head(Monocytes_mouse_2, n = 10)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Nog 9.502314e-11 0.3833897 0.333 0.020 1.297541e-06
## Nuf2 3.420788e-08 0.3746294 0.333 0.027 4.671086e-04
## Espl1 3.704791e-07 0.3686375 0.333 0.032 5.058893e-03
## Cd5 1.054815e-06 0.3659221 0.333 0.034 1.440349e-02
## Sfmbt2 1.310108e-06 0.3626703 0.333 0.034 1.788952e-02
## Crmp1 1.440750e-06 0.3653796 0.333 0.035 1.967344e-02
## Wbp1 2.018414e-06 0.4423256 0.444 0.061 2.756144e-02
## Zfp930 8.394765e-06 0.4700501 0.333 0.040 1.146305e-01
## Cbr4 2.050164e-05 0.4826788 0.556 0.108 2.799499e-01
## Camk4 5.347287e-05 0.4656349 0.556 0.116 7.301720e-01
write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Monocytes_mouse_2.csv")
Macrophage_cell_deseq <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(Macrophage_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophage_cell_deseq.csv")
saline_cellchat <- createCellChat(object = saline_mouse_1_analyzed, group.by = "ref.data")
## The `data` slot in the default assay is used. The default assay is SCT
## The `meta.data` slot in the Seurat object is used as cell meta information
## The cell groups used for CellChat analysis are B cells DC Endothelial cells Epithelial cells Fibroblasts ILC Macrophages Monocytes Neutrophils Stromal cells
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)
CellChatDB.use <- CellChatDB
saline_cellchat@DB <- CellChatDB.use
saline_cellchat <- subsetData(saline_cellchat)
saline_cellchat <- identifyOverExpressedGenes(saline_cellchat)
saline_cellchat <- identifyOverExpressedInteractions(saline_cellchat)
saline_cellchat <- computeCommunProb(saline_cellchat)
saline_cellchat <- filterCommunication(saline_cellchat, min.cells = 10)
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells: ILC Neutrophils
saline_cellchat <- computeCommunProbPathway(saline_cellchat)
saline_cellchat <- aggregateNet(saline_cellchat)
groupSize <- as.numeric(table(saline_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(saline_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(saline_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
# See interaction of each cell type with the another
mat1 <- saline_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
mat2[i, ] <- mat1[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}
##Visualize signalling pathways in each cell
# Identify signaling pathways showing significant communications
saline_cellchat@netP$pathways
## [1] "COLLAGEN" "APP" "MIF" "THBS" "VEGF"
## [6] "FN1" "LAMININ" "MHC-II" "ICAM" "ITGAL-ITGB2"
## [11] "COMPLEMENT" "PECAM1" "CD52" "CD22" "CD45"
## [16] "JAM" "PARs" "CHEMERIN" "SEMA4" "CXCL"
## [21] "EPHA" "PDGF" "GAS" "HSPG" "FGF"
## [26] "SEMA3" "AGRN" "PTPRM" "CDH5" "SEMA7"
## [31] "NOTCH" "SEMA6" "IL16"
# Evaluate one of the pathways
pathways.show <- c("NOTCH")
vertex.receiver = seq(2,5) # a numeric vector.
?netVisual_aggregate(saline_cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(saline_cellchat, signaling = pathways.show, layout = "chord")
# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(saline_cellchat, signaling = pathways.show, color.heatmap = "Reds")
netAnalysis_contribution(saline_cellchat, signaling = pathways.show)
pairLR.MHC <- extractEnrichedLR(saline_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair
# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(saline_cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)
## [[1]]
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(saline_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)
# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(saline_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)
saline_cellchat <- netAnalysis_computeCentrality(saline_cellchat, slot.name = "netP")
# the slot 'netP' means the inferred intercellular communication network of signaling pathways
# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(saline_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)
ht1 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "incoming")
ht1 + ht2
selectK(saline_cellchat, pattern = "outgoing")
nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "outgoing", k = nPatterns)
# river plot
netAnalysis_river(saline_cellchat, pattern = "outgoing")
# dot plot
netAnalysis_dot(saline_cellchat, pattern = "outgoing")
selectK(saline_cellchat, pattern = "incoming")
nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "incoming", k = nPatterns)
# river plot
netAnalysis_river(saline_cellchat, pattern = "incoming")
# dot plot
netAnalysis_dot(saline_cellchat, pattern = "incoming")
saline_cellchat <- computeNetSimilarity(saline_cellchat, type = "functional")
saline_cellchat <- netEmbedding(saline_cellchat, type = "functional")
## Manifold learning of the signaling networks for a single dataset
#> Manifold learning of the signaling networks for a single dataset
saline_cellchat <- netClustering(saline_cellchat, type = "functional")
## Classification learning of the signaling networks for a single dataset
#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(saline_cellchat, type = "functional", label.size = 3.5)
organism = org.Mm.eg.db
# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophage_cell_deseq.csv", header=TRUE)
# we want the log2 fold change
original_gene_list <- df$avg_log2FC
# name the vector
names(original_gene_list) <- df$X
# omit any NA values
gene_list<-na.omit(original_gene_list)
# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)
#gene set enrichment function for B cell cluster
gse <- gseGO(geneList=gene_list,
ont ="ALL",
keyType = "SYMBOL",
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = organism,
pAdjustMethod = "none",
eps = 0.0)
require(DOSE)
gse_dataframe <- as.data.frame(gse)
## count the gene number
gene_count<- gse@result %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)
## merge with the original dataframe
gse_dataframe<- left_join(gse@result, gene_count, by = "ID") %>% mutate(GeneRatio = count/setSize)
write.csv(gse_dataframe, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm+bcg/gse_dataframe.csv")
ggplot(gse_dataframe[c(3,9,17,21,27,29,41,58,59,63,89,102,109,118,136,139,145,163),],
aes(x = GeneRatio, y = Description)) +
geom_point(aes(size = GeneRatio, color = p.adjust)) +
theme_bw(base_size = 14) +
scale_colour_gradient(limits=c(0, 0.10), low="red") +
ylab(NULL) +
ggtitle("GO pathway enrichment")